using LinearAlgebra, Plots, Printf, SparseArrays
This implementation is again in Julia. Since tridiagonal matrices are so prevalent and important, Julia has a special data type for them. If you are using python or Matlab you want to use the spdiags command. For example, in python:
import numpy as np
from scipy.sparse import spdiags
m = 4
data = np.array([-np.ones(m), np.zeros(m) , np.ones(m)]);
diags = np.array([-1, 0, 1])
A = spdiags(data, diags, m, m)
A.toarray() # just to see what it looks like
h = 0.01
m = convert(Int64,1/h)-1;
k = 0.1
T = 10.
S = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) #increase dim by 1
a = 1.0;
S |> Array;
For $a > 0$
$$ \begin{cases} u_t + a u_{x} = 0,\\ u(x,0) = \eta(x),\\ u(0,t) = u(1,t). \end{cases} $$# Initial condition
η = x -> exp.(-20*(x .-1/2).^2)
# Initial condition chosen so that u(x,t) = sin(2*pi*(x - t)), if a = 1
# η = x -> sin.(2*pi*x)
# u = (x,t) -> sin.(2*pi*(x.-t))
#1 (generic function with 1 method)
Build MOL matrix
S = sparse(S) # Need to convert S to a new data type to allow new entries
S[1,end] = -1
S[end,1] = 1
S |> Array
100×100 Matrix{Float64}:
0.0 1.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 -1.0
-1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 -1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 -1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 -1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 -1.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
⋮ ⋮ ⋱ ⋮
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 … -1.0 0.0 1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 1.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 1.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 1.0
1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0
k = h^2
B = I - (a*k/(2h)*S)
100×100 SparseMatrixCSC{Float64, Int64} with 300 stored entries:
⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈
⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀
⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
# If solution is known
# plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
U = B*U
if mod(i-1,tb) ≈ 0.0
IJulia.clear_output(true)
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
# If solution is known
# plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
end
end
plot()
anim = Animation()
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)
frame(anim)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
U = B*U
if mod(i-1,tb) ≈ 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)
frame(anim)
end
end
gif(anim,"advection_periodic.gif")
┌ Info: Saved animation to │ fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/advection_periodic.gif └ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114
# Initial condition
η = x -> exp.(-20*(x .-1/2).^2)
# Initial condition chosen so that u(x,t) = sin(2*pi*(x - t)), if a = 1
#η = x -> sin.(2*pi*x)
#u = (x,t) -> sin.(2*pi*(x.-t))
#3 (generic function with 1 method)
h = 0.01
m = convert(Int64,1/h)-1;
T = 10.
S = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) #increase dim by 1
a = 1.0;
S = sparse(S) # Need to convert A₀ to a new data type to allow new entries
S[1,end] = -1
S[end,1] = 1
C₀ = sparse(Tridiagonal(fill(0.5,m),fill(0.0,m+1),fill(0.5,m)))
C₀[1,end] = .5
C₀[end,1] = .5
k = h # need this to avoid too much numerical dissipation
B = sparse(C₀ - (a*k/(2h))*S);
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
# If solution is known
# plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
U = B*U
if mod(i-1,tb) ≈ 0.0
IJulia.clear_output(true)
plot(x, U, xaxis = [0,1], yaxis = [0,1],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
# If solution is known
# plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
end
end
plot()
anim = Animation()
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)
frame(anim)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
U = B*U
if mod(i-1,tb) ≈ 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[U[end],U[end]], label = "BCs", seriestype = :scatter)
frame(anim)
end
end
gif(anim,"advection_periodic_LF.gif")
┌ Info: Saved animation to │ fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/advection_periodic_LF.gif └ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114
Consider the method of lines discretization:
$$ U'(t) = - \frac{a}{2h} A U(t) + \frac{a}{2h} g(t)$$where
$$A = \begin{bmatrix} 0 & 1 \\ -1 & 0 & 1 \\ & -1 & 0 & 1 \\ & & \ddots & & \ddots\\ &&& -1 & 0 & 1 \\ &&& 1 & -4 & 3 \end{bmatrix}, \quad g(t) = \begin{bmatrix} g_0(t) \\ 0 \\ \vdots \\ 0 \end{bmatrix}.$$Forward Euler produces:
$$U^{n+1} = U^n - \frac{ak}{2h} A U^n + \frac{ak}{2h} g(t). $$We introduce a Lax-Friedrichs stabilization via a matrix $C$
$$ C =\begin{bmatrix} 0 & 1/2 \\ 1/2 & 0 & 1/2 \\ & 1/2 & 0 & 1/2 \\ & & \ddots & & \ddots\\ &&& 1/2 & 0 & 1/2 \\ &&& & 0 & 1 \end{bmatrix}.$$Using the approximation $U^n \approx C U^n + g(t)/2$ we find the method
$$U^{n+1} = C U^n - \frac{ak}{2h} A U^n + \left( \frac 1 2 + \frac{ak}{2h}\right) g(t). $$
# Initial and boundary conditions
η = x -> exp.(-20*(x .-1/2).^2)
g0 = t -> sin(4*t)+1.
#7 (generic function with 1 method)
η(0.0) |> display
g0(0)
0.006737946999085467
1.0
h = 0.01
a = 1.0;
m = convert(Int64,1/h)-1;
k = h/(a)
T = 10.
A = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)); #increase dim by 1
A = sparse(A)
A[end,end-2:end] = [1,-4,3];
C = Tridiagonal(fill(0.5,m),fill(0.0,m+1),fill(0.5,m))
C = sparse(C)
C[end,:] *= 0.0
C[end,end] = 1.0
B = sparse(C - (a*k/(2h))*A);
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
U = B*U
U[1] += (1/2 + a*k/(2h))*g0(t-k)
if mod(i-1,tb) ≈ 0.0
IJulia.clear_output(true)
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
end
end
plot()
anim = Animation()
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter)
frame(anim)
fr = 100 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
U = B*U
U[1] += (1 + a*k/h)*g0(t)/2
if mod(i-1,tb) ≈ 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter)
frame(anim)
end
end
gif(anim,"advection_LF.gif")
┌ Info: Saved animation to │ fn = /Users/thomastrogdon/Dropbox (uwamath)/Teaching/586/amath-586-2022/notebooks/advection_LF.gif └ @ Plots /Users/thomastrogdon/.julia/packages/Plots/zozYv/src/animation.jl:114
# Initial and boundary conditions
η = x -> exp.(-40*(x .-1/2).^2)
g0 = t -> sin(4*t) + 1
# Initial condition chosen so that u(x,t) = sin(2*pi*(x - t)), if a = 1
# η = x -> sin.(2*pi*x)
# g0(t) = sin.(-2*pi*t)
# u = (x,t) -> sin.(2*pi*(x.-t))
#11 (generic function with 1 method)
h = 0.01
a = 1.0;
m = convert(Int64,1/h)-1;
k = h/(a)
T = 10.
A = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) |> sparse;
D = Tridiagonal(fill(1.0,m),fill(-2.0,m+1),fill(1.0,m)) |> sparse;
vec = [1,-4,3]
A[end,end-2:end] = vec; # same as Lax-Friedrichs
# Construct a backward second-order approximation of the second derivative
D[end,:] *= 0.0
D[end,end-3:end-1] += (vec[1]/2)*[-.5,0,.5]
D[end,end-2:end] += (vec[2]/2)*[-.5,0,.5]
D[end,end-2:end] += vec[3]*vec/4;
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
# If solution is known
# plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter)
fr = 1000 #frames/unit time
tb = convert(Int64,ceil(n/(fr*T)))
for i = 2:n+1
t += k
U = U - (a*k/(2h))*(A*U) + (a^2*k^2/(2h^2))*(D*U)
U[1] += (a*k/(2h) + (a^2*k^2/(2h^2)))*g0(t-k)
if mod(i-1,tb) ≈ 0.0
IJulia.clear_output(true)
plot(x, U, xaxis = [0,1], yaxis = [-1,2],lw=3,label = @sprintf("u(x,t), t = %1.2f",t))
# If solution is known
# plot!(x, u(x,t), xaxis = [0,1], yaxis = [-1,2],lw=1,label = @sprintf("u(x,t), t = %1.2f",t))
plot!([0,1],[g0(t),U[end]], label = "BCs", seriestype = :scatter) |> IJulia.display
end
end
# Initial and boundary conditions chosen so that u(x,t) is the solution
F = x -> (exp.(-(x .- 1).^2) .+ 1).*sin.(2*pi*x)
a = 1.0
η = x -> F(x)
g0 = t -> F(-a*t)
u = (x,t) -> F(x.-a*t)
#19 (generic function with 1 method)
h = 0.2
p = 7
out = fill(0.0,p)
for i = 1:p
h = h/2
a = 1.0;
m = convert(Int64,1/h)-1;
k = h/(a)
T = 10.
A = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) |> sparse;
D = Tridiagonal(fill(1.0,m),fill(-2.0,m+1),fill(1.0,m)) |> sparse;
vec = [1,-4,3]
A[end,end-2:end] = vec; # same as Lax-Friedrichs
# Construct a backward second-order approximation of the second derivative
D[end,:] *= 0.0
D[end,end-3:end-1] += (vec[1]/2)*[-.5,0,.5]
D[end,end-2:end] += (vec[2]/2)*[-.5,0,.5]
D[end,end-2:end] += vec[3]*vec/4;
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
for i = 2:n+1
t += k
U += -(a*k/(2h))*(A*U) + (a^2*k^2/(2h^2))*(D*U)
U[1] += (a*k/(2h) + (a^2*k^2/(2h^2)))*g0(t-k)
end
out[i] = maximum(abs.(U - u(x,T)))
end
out[1:end-1]./out[2:end]
6-element Vector{Float64}:
7.382313823819206
7.841350172193508
7.960061016385922
7.989999195538282
7.997501306853115
7.999007598621097
η = x -> x.*exp.(-10*(x .- 1).^2)
η(0.00001)
η = x -> exp.(-20*(x .-1/2).^2)
g0 = t -> sin(4*t)^2.
h = 0.2
h = h/2^10
a = 1.0;
m = convert(Int64,1/h)-1;
k = h/(a)
T = 10.
A = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) |> sparse;
D = Tridiagonal(fill(1.0,m),fill(-2.0,m+1),fill(1.0,m)) |> sparse;
vec = [1,-4,3]
A[end,end-2:end] = vec; # same as Lax-Friedrichs
# Construct a backward second-order approximation of the second derivative
D[end,:] *= 0.0
D[end,end-3:end-1] += (vec[1]/2)*[-.5,0,.5]
D[end,end-2:end] += (vec[2]/2)*[-.5,0,.5]
D[end,end-2:end] += vec[3]*vec/4;
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
for i = 2:n+1
t += k
U += -(a*k/(2h))*(A*U) + (a^2*k^2/(2h^2))*(D*U)
U[1] += (a*k/(2h) + (a^2*k^2/(2h^2)))*g0(t-k)
end
U_true = U;
h = 0.2
p = 7
out = fill(0.0,p)
for i = 1:p
h = h/2
a = 1.0;
m = convert(Int64,1/h)-1;
k = h/(a)
T = 10.
A = Tridiagonal(fill(-1.0,m),fill(0.0,m+1),fill(1.0,m)) |> sparse;
D = Tridiagonal(fill(1.0,m),fill(-2.0,m+1),fill(1.0,m)) |> sparse;
vec = [1,-4,3]
A[end,end-2:end] = vec; # same as Lax-Friedrichs
# Construct a backward second-order approximation of the second derivative
D[end,:] *= 0.0
D[end,end-3:end-1] += (vec[1]/2)*[-.5,0,.5]
D[end,end-2:end] += (vec[2]/2)*[-.5,0,.5]
D[end,end-2:end] += vec[3]*vec/4;
n = convert(Int64,ceil(T/k))
x = h:h:1 # include right end point
U = η(x)
t = 0.0
for i = 2:n+1
t += k
U += -(a*k/(2h))*(A*U) + (a^2*k^2/(2h^2))*(D*U)
U[1] += (a*k/(2h) + (a^2*k^2/(2h^2)))*g0(t-k)
end
U_true_p = U_true[end:-2^(10-i):1]
out[i] = maximum(abs.(U[end:-1:1] - U_true_p))
end
out[1:end-1]./out[2:end]